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Monte Carlo simulation of joint density of states of two continuous spin models using 

Wang-Landau- Transition-Matrix Algorithm 

Shyamal Bhaio and Soumen Kumar RojO 
Department of Physics, Jadavpur University, Kolkata-700032, India. 

Monte Carlo simulation has been performed in one-dimensional Lebwohl-Lasher model and two 
dimensional XY-model using the Wang-Landau and the Wang-Landau- Transition-Matrix Monte 
Carlo methods. Random walk has been performed in the two-dimensional space comprising of 
energy-order parameter and energy-correlation function and the joint density of states (JDOS) 
were obtained. From the JDOS the order parameter, susceptibility and correlation function are 
calculated. Agreement between the results obtained from the two algorithms is very good. 

PACS: 64.60.De; 61.30.-v; 05.10.Ln 

Keywords: Wang Landau, Transition Matrix. Joint Density of states 



(N 

o 

B 

I 

-(— > 

-(— > 



X3 
O 

o 



> 

00 

o 
o 



X 



I. INTRODUCTION 



' During the last couple of years or so a number of Monte- 
, Carlo (MC) algorithms have been proposed which di- 
' rectly determine the density of states (DOS) of a sys- 
' tem. One of these is the Transition Matrix Monte Carlo 
(TMMC) algorithm developed by Oliveria et al^ and sub- 
I sequently generalized by Wang and co-workers^. In this 
I algorithm, during a random walk in the energy space, one 
' keeps a record of the transitions between the microstates 
of the system. The entire history of the transitions is 
' then used to obtain the density of states of the system. 
I More recently Wang and Landau'^ proposed another al- 
gorithm which goes by their name (WL) and has drawn 
^ wide attention of investigators. This algorithm employs 
. a method of flat histograms, while a random walk is per- 
' formed in the energy space and estimates the DOS of the 
system by using an iterative scheme. Both methods, the 
TMMC and WL, depend on broad sampling of the phase 
space and are easy to implement. A knowledge of the 
DOS of a system as a function of energy, il{E) enables 
one to calculate the partition function Z by a simple 
Boltzmann reweighting: Z(j3) — 



^^n{E)e-^'^, where 



E 

P is the inverse temperature 1/T. With a knowledge of 
the partition function one can calculate the averages of 
thermodynamic quantities which are directly related to 
energy. It has been established that while the TMMC 
method gives more accurate estimation of the DOS, the 
WL algorithm is more efficient in sampling the phase 
space. 

Shell et al proposed^ an algorithm which is an amalga- 
mation of the two algorithms and utilizes the benefits of 
each. This algorithm, now known as the Wang Landau 
Transition Matrix (WLTM) Monte Carlo algorithm, is at 
the same time efficient and accurate. Shell et al applied 
the algorithm to a two-dimensional Ising model and a 
Lennard- Jones fluid. More recently Ghulghazarya et a& 
have applied this method to simulate protein and pep- 



tide. 

The sampling of the phase space which is done in each 
of the above methods need not be restricted to the eval- 
uation of the density of states as a function of energy 
alone. One can determine the DOS or to be more spe- 
cific, the Joint Density of States (JDOS) with a sub- 
stantially more book-keeping. The JDOS ^{E, (p) is a 
function of some variable (which can be an order pa- 
rameter, spin-spin correlation or any other observable) 
besides the energy E. This partition function is deter- 
mined from: Z(^) = X! X! ^(^' <^)e~'^'^- The ensemble 

E <t> 

average of any function of at an inverse temperature /3 
is then given by 



^^/(0)f](ii;,0)e-^^ 



</(0,/3)> = 



Y,Y.^{E,4>)e 



I3E 



(1) 



Most of the investigators have so far worked on the de- 
termination of DOS or JDOS in discrete systems using 
the three above mention algorithms^'"— . Even in this 
domain the amount of work reported on the determina- 
tion of JDOS is relatively small. In an earlier paper— 
we have reported on the working of WLTM in two con- 
tinuous lattice-spin models. While the discrete systems 
like Ising or Potts model can be handled in a straight 
forward manner, the investigation of continuous systems 
are more tedious. The range of energy (and another ob- 
servable in case of two-dimensional random walk) needs 
to be discretized and several parameters appear which 
are to be chosen properly for the determination of DOS 
or JDOS. In the present communication we report the 
determination of JDOS using the WLTM algorithm in 
two continuous lattice spin models. The quantities eval- 
uated other than energy and specific heat are the or- 
der parameter, susceptibility and correlation function. 
One of the models is the one-dimensional Lebwohl-Lasher 
model^^, which is exactly solvable!^ and therefore allows 



us to check the accuracy of the results of our simula- 
tion. The other system is the two-dimensional XY-model 
where exact solutions are not available and we have com- 
pared the results with those obtained from the JDOS 
determined using the WL algorithm. The aim of the 
present work is to test the feasibility of the determina- 
tion of JDOS in continuous models using the WLTM al- 
gorithm. This turned out to be a some what difficult task 
as a formidable amount of computer memory is necessary 
even for systems of moderate size. 



II. THE DIFFERENT ALGORITHMS 

A. The Wang Landau algorithm 

We outline below the method of determination of the 
JDOS, fl{E, (j)) using the WL algorithm. In a system 
where E and </> are continuous variables, discretization 
in required to label the macrostates of the system. The 
ranges of E and 4> are divided into a large number of 
bins of width de and d^ respectively. Let Ei and cjjj 
be the mean energy and mean order parameter (or 
correlation function) corresponding to the /*'* bin of 
energy and the J*'' bin of the order parameter (corre- 
lation function) respectively. Then ri(/, J) denotes the 
number of microstates of the system having energy Ej 
and order parameter (correlation function) 0j or simply 
the degeneracy of the macrostate {I, J). Since n{I,J) 
is a very large number, it is convenient to work with its 
natural logarithm and we use g{I, J) — lnD,(I, J). 

We perform two dimensional random walk in the {E — (p) 
space. Initially we do not have any knowledge oi g{I, J), 
and set g{I, J) = for all values of / and J. Also, an 
histogram count H{I, J) of the states visited during the 
random walk is maintained. The WL algorithm gener- 
ates the JDOS profile, which progressively approaches 
the actual density of states of the system. The algorithm 
starts with some microstate of the system and successive 
microstates are generated by rotating one spin at a time. 
Let (/, J) and {K, L) be the macrostates before and after 
rotating the spin and the corresponding microstates are 
(i,j) and {k,l), where i & I,j & J and k € K,l G L. The 
transition probability from state {i,j) to (fc, I) is given by 



and 



H{K,L) = H{K,L) + 1 



(3b) 



P{iJ^k,l) 



fn(i,j) 



(2) 



Thus the acceptance probability of the new state is in- 
versely proportional to the current density of states. 
When the new state is accepted the density of states 
g{K,L) and the histogram count H{K,L) of the state 
{K, L) are modified as 



and when the new state is not accepted the old density of 
state g{I, J) and histogram count H{I, J) are modified 
as 



g(/, J) = g(/, J) + Inf 



and H{I, J) = H{I, J) + 1 



(4a) 



(4b) 



g{K,L)^g{K,L) + lnf 



(3a) 



Here / is a modification factor whose initial value 
was chosen to be equal to e. When the histogram is 
sufBciently flat (say, 80%) i.e., histogram count H{I,J) 
of each bin (/, J) is at least 80% of the mean histogram 

Hm = (Ef ' Ef ' H{I, J)) /M , M being the total 
number of bins then one iteration is said to be complete. 
Then the histogram is reset to zero for all values of 
/ and J and the modification factor / is reduced in 
some prescribed manner (we use Inf — s> lnf/2). A 
fresh iteration is started with the modified value of Inf 
and the old values of g{I, jys which were calculated 
in the previous iteration. One continues iterations 
with the same procedure until the modification factor 
becomes sufficiently small (say, 10~*). The error intro- 
duced in the joint density of states has been predicted 
to be proportional to \/{lnf) as is apparent from 
the theoretical work of Zhou and Bhatt^i^. This has 
been tested for a number of discrete and continuous 
models and the prediction has been found to be correcti^. 



B. The Transition Matrix Monte Carlo algorithm 

The TMMC algorithm is an efficient algorithm in which 
one directly calculates the density of states and was first 
proposed by Oliveria et al. in the year 1996. If the 
transition probability from a microstate {i,j) to another 
microstate (fc, I) is t{i,j; k, I) and that from a macrostate 
(/, J) to a macrostate {K, L) is T{I, J; K, L) then 

with the following conditions: 

J2t{i,j;k,l) = l, t{i,j;k,l)>0 (6a) 

k,l 

and ^T(/, J;i^,i) = 1, T{I,J;K,L)>0 

K,L 

(6b) 
If T{K, L; I, J) be the reverse transition then one can 
write 



T{I, J; K, L) ^{K, L) 0=1 jej keK leL 
TiK,m,J)= mi, J) 5:^^^t(fc,/;z,j) 

keK IGL iGl j&J 



(7) 



The transition probability t{i,j] k, I) actually is a product 
of two probabilities, 



tihj] k, I) = a{i,j; k, l)P{i, j; k, I) 



(8) 



Here a(i,j;k,l) is the probability of the move (j,j) — > 
(fc, I) being proposed and depends on the type of the 
Monte Carlo moves and while P{i,j;k,l) is the proba- 
bility of proposed move being accepted and depends on 
the configurations (i,j) and {k,l). One can choose any 
value of P and an infinite temperature transition prob- 
ability (Too) can be chosen for which P{i,j; k,l) = 1 for 
all i,j and k,l. This is particularly easy to understand if 
one considers the Metropolis algorithm at infinite tem- 
perature. This process does not affect the JDOS to be 
determined since it is independent of temperature. Using 
equations ([7]) and ([5]) we can write: 

Too{I,J;K,L) _ fljK, L) iei je.J keK leL 

keK leL iei jeJ 

(9) 
Again for symmetric moves (single spin flip dynamics) 
a{i,j;k,l) = a{k,l;i,j) and the summation terms on 
right hand side of equation ^ drops out and the equa- 
tion can be simplified as 



T^{I,J;K,L) _n{K,L) 



T^{K,L;I,J) nil,. J) 



(10) 



This equation relates the JDOS with the infinite tem- 
perature transition probabilities. Thus from the knowl- 
edge of infinite temperature transition probabilities one 
can estimate the JDOS. Now, one needs to calculate 
Too(I,J;K,L) which can be done by keeping a record 
of moves in the form of a matrix, called C-matrix, 
C{I,J;K,L) for all proposals {I, J) — > {K,L), during 
the random walk. Initially, we set C{I, J; K,L) = for 
all I, J and K,L. At infinite temperature all the proposed 
moves are accepted so whenever a move is proposed we 
update the C-matrix as 



C(/, J; K, L) = C{I, J; K,L) + l 



(11) 



Once the construction of C-matrix is started we never 
reset it to zero, because the C-matrix keeps the detail 
history of the transitions in the system. The current es- 
timate of the infinite temperature transition probability 



[To. 



IS 



T^iI,J;K,L) 



C{I,J;K,L) 



(12) 



K 



Where the sum extends over all K and L and the tilde 
indicates the estimate. From equation ((T0| and (J12p one 
can determine the joint density of states. But the equa- 
tion (jlOp is an over specified problem since for a sys- 
tem with N macrostate having N unknown quantities 
n{I, J) there are N{N — l)/2 such equations. To calcu- 
late n(I, J) one needs to minimize the total variance 



2 

^tot 



with 



E 

I,J;K,L 



[S{I,J)-S{K,L) + ln(^^^ 



{I,J;K,L) 
(K,L-I,J) 



2 
'^IJKL 



(13) 



(tIjkl^C{I,J;K,L)-^+H{I,J) 



'^+C{K, L; I, J)-^+H{K, L)-^ 
(14) 



Here we have used S{I, J) = lnVt{I, J) and H{I, J) = 
'Y^j^ ^C{I,J;K,L). In order to ensure that the ran- 
dom walker visits all macrostates in the region of interest 
one considers a uniform ensemble, where all macrostates 
are equally probable. The probability of occurrence of 
a given microstate («, j) is therefore proportional to the 
multiplicity of the macrostate (/, J) where i G / and 
i e J. So the probability of acceptance of a move 
{i,j) -^ {k,l) is given by 



p{i,j — >■ k,l) — min 1 



»(/, J) 

n{K, L) 



(15) 



Since the JDOS are not known a priori, the acceptance 
criteria can be written as 



p{i,j -^ k,l) — min 1 



To^{K,L;I,J) 
f^{I,J;K,L) 



(16) 



where we have used equation (jlOp So using the above 
acceptance probability together with equation (llOp and 
minimizing the total variance as discussed above, one can 
generate the profile of the joint density of states. The 
TMMC algorithm gives accurate value of JDOS since it 
stores and uses the entire history of the transitions dur- 
ing the random walk. But this method has a drawback 
that the convergence is not guaranteed and a significant 
amount of CPU time is necessary. On the other hand the 
WL method ensures that all the macrostates are visited 
rather efficiently and more quickly. As proposed by Shell 
and co-workers- in the WLTM method we combine the 
TMMC and WL algorithms in order to get the benefit of 
both the algorithms. This is discusses in the next section. 



C. Combination of Wang-Landau and Transition 

Matrix Monte Carlo algorithm: The WLTM 

algorithm 

In the WLTM algorithm, one efficiently combines the 
WL and TM Monte Carlo algorithms. In this algorithm, 
the phase space is sampled via the acceptance criteria 
given by equation ([15]) as in the WL method. Also 
g{I, J) and H{I, J) are updated in the same fashion as 
in the WL algorithm and in addition of these, a record 
of transitions between macrostates is kept in a C-matrix, 
which is never zeroed during the simulation. At the end 
of the simulation, from the knowledge of the C-matrix 
and using equation (jlOp the joint density of states is 
determined by minimizing the total variance. 



III. THE MODELS USED IN THE PRESENT 
WORK 



then — ^(a) have equal contribution as that of 77(a) to- 
wards the order parameter. So a vector order parameter 
is inadequate for the system and a symmetric traceless 
tensor is used as the order parameter. This tensor is cho- 
sen in such a way that it is zero in the high temperature 
isotropic phase and is unity in the fully ordered phase. 
The order parameter tensor is given as 



1 



N 



i=l 



h^ 



(18) 



where nf is the z*'* component of the unit vector n, which 
points along the spin at position x^a. N is the number of 
molecules in the system. The order S ^< q- > prevailing 
in the system can be written as 

^-<.^>-^<^T.Q^-l> (19) 

The second rank spin-spin correlation function p{r) is 
defined as 



We have determined the JDOS and other thermody- 
namic quantities such as energy, specific heat, order 
parameter, correlation function etc. performing two 
dimensional random walk in E — (p space using WL 
and WLTM algorithms for two continuous lattice spin 
models. We give below a brief description of the models 
and the related quantities we have measured. 



A. The 1-d Leb^vohl-Lasher model 

This model is a linear array of three dimensional spins 
{d — l,n — S) interacting with nearest neighbors via the 
a potential 



p{r) ==< P2[cos0[r)) > 



(20) 



Vij = -P2{cos9ij) 



(17) 



Where P2 is the second Legendre polynomial and 6ij is 
the angle between two nearest neighbor spins i and j. 
Spins are headless, i.e. it has 0(3) as well as Z2 sym- 
metry. This model represents one dimensional nematic 
liquid crystal and does not exhibit any finite temperature 
order disorder phase transition. 

The order parameter is a quantity which describes the 
amount of order prevailing in a system. Since in a ne- 
matic phase the ordering is the orientational, the order 
parameter quantifies the amount of orientational order 
present in the system. In a nematic liquid crystal the 
molecules are on the average aligned along a particular 
direction n called the director. These molecules in gen- 
eral have equal probability of pointing parallel and anti- 
parallel to any given direction. If ^(a) is the direction of 
a molecular axis of a molecule situated at position x = a 



Where 9{r) is the angle between two spins separated 
by r lattice spacing r. The order parameter S and 
limr_>oo p{r) should vanish in the thermodynamic limit. 
However due to the finite size effect in systems of finite 
size both quantities have some small non zero value. The 
exact expression for the correlation function is given byi^ 



PN{r) 



3^-1/2^-1(^1/2) _ 3^-1 _i 

4 ^ ^ 4 2 



(21) 



where K — 3/2T. D is Dawson function—, 



D{x) = exp[~x ] / du exp[u 
Jo 



B. The 2-d XY model 

In this model planar spins placed at the sites of a pla- 
nar square lattice interact with nearest neighbours via a 
potential. 



F(%) = 2{l-[cos2(0,,/2)]} 



(22) 



where 9ij is the angle between the nearest neighbours i, j . 
[This particular form of the interaction, rather than the 
more conventional ~cos{9ij) form, was chosen by Do- 
many et. ai^ to enable them to modify the shape of the 
potential easily, which led to what is now known as the 
modified XY- model ]. The XY-model is known to ex- 
hibit a quasi-long-range-order disorder transition which 
is mediated by unbinding of topological defects as has 
been described in the seminal work of Kosterlitz and 



Thoule a ^'^^. The XY-model has also been the subject of 
extensive MC simulation over last few decades and some 
of the recent results may be found in^°. 
In this model, the orientational order parameter is de- 
fined as follows: let n be the unit vector (called the di- 
rector) in the direction of maximum order prevailing in 
the system and s be the spin vector of unit magnitude 
then the order parameter S is given by 



S =< n • s > =< cos 6 > 



(23) 



where (p is the angle between the director and the spin. 
The spin- spin correlation function for the 2d-XY model 
is defined as 



p{r) 



=< cos Oir) > 



(24) 



Where 9{r) is the angle between two spins separated by 
r lattice spacing. 



IV. COMPUTATIONAL DETAILS: 

Two dimensional random walk has been performed in 
E — S and E — p space for one dimensional Lebwohl- 
Lasher model and two dimensional XY model. Since both 
of these models are continuous, we discretize the system 
by dividing the whole energy range as well as the order 
parameter (or correlation function) range into a number 
of bins having widths de and d^ respectively. We simu- 
lated 2d XY system for lattice sizes 5, 10, 15 and 20 and 
Id LL system for lattice sizes 80, 160 and 220. The 2d 
XY system can have energy between and 2N'^ but we 
simulated the system for the energy range 3.0 to 2N'^ to 
avoid trapping of the random walker in low energy states 
as these are scarcely visited during the simulation. Simi- 
larly the Id LL system can have energy between —N and 
N/2 but we simulated the system for the energy range 
—N, with a small energy cut near the ground state. 
We have chosen de to be 0.1 for the Id LL model and 0.2 
for the XY model, d^ was chosen to be 0.01 for both the 
models while performing random walk in E — S space as 
well as in the E — p space. In the Id LL model, each 
spin have three components (h with i — 1,2,3) specify- 
ing the three direction cosines of the spin. We have taken 
a random initial configuration and a new configuration is 
generated by rotating any one of the spins randomly us- 
ing the prescription k = li + p * ri (for i=l,2,3), where 
fj is a random number between -1 and 1. In the 2d XY 
model, each spin is specified by two direction cosines and 
new configurations are generated by rotating the spin in 
the same manner as above. Initially, we do not have any 
prior knowledge about the DOS, so we set g{I, J) = for 
all values of the macrostates / and J. A new microstate 
(fc, I) is generated from the old one (i, j) by rotating one 
spin at a time and the acceptance probability p{i,j; k, I) 



is given by equation ^. When the new state lies in the 
same macrostate as the old one the state is accepted. 
A histogram count is recorded in an array H{I, J) and 
whenever a state is proposed we update the C-matrix as 
C(/, J; K, L) = C(/, J;K,L) + 1. The random walk is 
continued till the histogram becomes flat (say 80%). We 
point out that in the case of 2-d random walk all the bins 
are not visited uniformly as this would need an enormous 
number of Monte Carlo Sweeps (MCS). We first sam- 
pled the system for about 10® MCS which is called the 
'pre-production run'. During the pre-production run we 
marked by '1' the bins which are visited at least 80% of its 
average value and by '0' those which are visited less than 
80% of the average value or are not visited at all. During 
the 'production run' we checked the histogram flatness 
only for those bins which were marked '1' . The average 
histogram is calculated by considering only those bins 
which are visited at least once discarding the bins which 
are not visited at all. Some bins marked '0' may qualify 
for the mark '1' during the production runs. When the 
histogram gets flat the modification factor is changed as 
Inf -^ Inf /2 and the histogram count is reset to zero 
but the C-matrix is never reset to zero. The iteration 
is continued with the new modification factor and the 
same procedure is repeated until the modification factor 
becomes as small as 10~^ for LL model and 10~^ for XY 
model. At the end of the simulation the JDOS is calcu- 
lated from the knowledge of the C-matrix by minimizing 
the variance in equation (1131) . This variance is minimized 
using the Dowhill-Simplex method^i. 
It may be noted that for the continuous system there are 
a large number of macrostates / and J. As a result, the 
number of elements in the four dimensional C-matrix is 
enormously large. To simulate the system a huge amount 
of computer storage (RAM) is required which is beyond 
of our computer resources. But most of the elements of C- 
matrix are zero since a transition from a given macrostate 
(/, J) to all macrostates {K, L) is not possible. In gen- 
eral it is found that for the transition (J, J) — >■ {K, L), the 
possible nonzero values of K and L lie within the range 
I — III to I+ni and J~n2 to J+n2 respectively, where ni 
and 712 are integers. Some extra conditions are imposed 
for the low values of / and J. In our case the values of ni 
and 71.2 are 5 and 1 respectively. So to avoid the problem 
with storage we transformed K to K' — K — I + ni and 
L to L' = L — J + n2 and whenever the original values 
of K and L are necessary we make the reverse transfor- 
mations. 

The JDOS is obtained by minimizing the variance using 
the Downhill Simplex method^i. The variance is a func- 
tion of more than one independent variable; in fact it is 
a function of a large number of independent variables. It 
is also very difhcult to minimize the multivariable func- 
tion of such a large number of unknown variables using 
this method. For example, to minimize a function of N 
variables (which constitutes a point in an N-dimensional 



space) it requires A^ + 1 points in that space. One of 
these points is chosen as the starting point Pq and other 
points are obtained by Pi = Pa + Ae^, where A is a con- 
stant and is taken to be 0.1% of the value of Pq, and 
the Ci s are N unit vectors. In this method we need a 
two dimensional matrix (p) of extents N + 1 and A'^. For 
a continuous model, the values of N is very large con- 
sequently the p matrix becomes very large. We had to 
divide the entire two dimensional surface g{I, J) into a 
large number of smaller segments. The energy-order pa- 
rameter surface g{I, J) of LL model was divided into 79, 
159, 218 segments for the lattice sizes 80, 160 and 220 re- 
spectively. Each segment consists of 100 order parameter 
bins and 10 energy bins. Similarly, the energy-correlation 
function surface of LL model is divided into the same 
number of segments as before each having 140 correla- 
tion function bins and 10 energy bins. We divided the 
energy-order parameter surface of XY model into 224, 
984, 2132 and 3980 segments for lattice sizes 5, 10, 15 
and 20 respectively. Each segment contains 100 order 
parameter bins and 14, 12, 12, 10 energy bins for lattice 
sizes 5, 10, 15 and 20 respectively. Similarly the energy 
correlation function surface of the XY model was also 
divided into the same number of segments. Each of the 
segments contains 100 correlation function bins and 10 
energy bins. Minimization of the variance given by equa- 
tion (fT3|) was carried out separately over each segment 
using the Downhill Simplex method and the JDOS for 
the entire surface is obtained by connecting the JDOS of 
each segment obtained by the above mentioned minimiz- 
ing method. With the above procedure we were able to 
calculate the JDOS of energy-order parameter or energy- 
correlation function by minimizing the variance as stated. 
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FIG. 2: Plot of the logarithm of density of states, g(E) against 
system energy per particle for 1-d LL model of lattice size L — 
80, 160 and 220 obtained using WL and WLTM algorithms. 
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FIG. 1: The logarithm of JDOS for the L = 80 1-d LL model 
obtained from the WLTM method is plotted against the en- 
ergy and order parameter. 



FIG. 3: The variation of orientational order parameter (S) 
fo the 1-d LL model for different lattice sizes is plotted with 
temperature, obtained from both WL and WLTM algorithms. 



V. RESULTS AND DISCUSSION 

A. The 1-d LL model: 

In this model simulation was carried out for systems 
having linear dimension L = 80, 160 and 220. The 
JDOS obtained by using the WLTM algorithm is shown 
in figure 1 as a function of energy and order parameter 
for L = 80. In figure 2 we have plotted the DOS 
(actually its logarithm, g{E)) for the three lattices as a 
function of energy per particle obtained from the JDOS 
of figure 1 by using both WL and WLTM methods. 
From the knowledge of g{E) one can calculate the 
average energy, specific heat etc. To determine the order 
parameter 5'(T), one needs to perform 2d random walk 
in energy-order parameter space. Figure 3 shows the 
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FIG. 4: The Susceptibility is shown as a function of temper- 
ature for the 1-d LL model of different lattice sizes. Data 
obtained from WL algorithm as well as from WLTM algo- 
rithm are plotted in this graph. 



FIG. 6: Correlation function as a function of lattice spacing 
r for seven different temperatures. Both WL and WLTM 
algorithms have been used. Comparison has been made with 
the exact results (the dots) of-^ obtained using equation 1211 
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FIG. 5: Variation of correlation function p(r, T) with tem- 
perature for lattice size L = 80 of 1-d LL model. 15 different 
values of r taken 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28 
and 30. The top most curve is for r = 2 and the lower curves 
are for other values of r given in the sequence above and in 
ascending order of r. The curves are the results obtained 
from the JDOS obtained using WL and WLTM algorithms. 
Comparison has been made with the exact results (the dots) 
obtained fronji^ using equation [2l1 for r = 2, 4, 6, 8, 10, 12, 
14 and 16. 



order parameter as a function of temperature for the Id 
LL model for different lattice sizes obtained using WL 
and WLTM algorithms. From the figure it is clear that 
S decreases rapidly with temperature and with increase 
in system size. This is expected since the Id LL model 
does not possess a true long range order at any finite 
temperature and the non-zero values of S obtained is due 
to the finite size effect. Consequently the Id LL model 
does not exhibit any finite temperature order-disorder 



FIG. 7: Logarithm of density of states vs energy per particle 
for the 2-d XY model of sizes 5x5, 10 x 10, 15 x 15 and 
20 x 20. The results obtained from both WL and WLTM 
algorithm are plotted. 



phase transition. In figure 4, susceptibility is plotted 
against temperature for different system sizes. The 
peaks become sharp with the increase of system size. 
This also supports the absence of order-disorder phase 
transition. 

The correlation function p(r^ T) has been plotted against 
temperature with r, the spacing between spins, in 
figure 5 for lattice size L=80 and 15 different values 
of r ranging from 2 to 30. All data are obtained from 
2-d random walk in E — p space using both WL and 
WLTM algorithms. It may be noted that, we had to 
run one simulation for each value of r. With increase 
in temperature, the correlation between spins for a 
given value of r decreases. Again p(r, T) is plotted 
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FIG. 8: The orientational order parameter is shown against 
temperature T for 2-d XY model of different lattice sizes using 
both WL and WLTM algorithms. 
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FIG. 9: The susceptibility of 2-d XY model for four liner 
lattice sizes L = 5, 10, 15 and 20 against temperature is 
plotted. All the data are from the results of two dimensional 
random walk in _B — S* space using WL and WLTM algorithm. 
The height of the peak increases with the increase in system 
size. It is also noted that the position of peak shifts towards 
lower temperature with the increase of system size. 



FIG. 10: Variation of correlation function p(r, T) with tem- 
perature for lattice size L = 20 of the 2-d XY model. 10 
different values of r taken 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10. The 
top most curve is for r — 1 and the lower curves are for other 
values of r given in the sequence above and in ascending val- 
ues of r. The curves are the results we obtained from joint 
density of states obtained using the both WL and WLTM 
algorithms. 
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FIG. 11: Correlation function as a function of lattice spacing 
r with eight different temperatures as parameter is plotted. 
Both the WL and WLTM algorithms were used. 



The 2-d XY model 



as a function of lattice spacing r for seven different 
temperatures in figure 6. p{r,T) decreases very rapidly 
with lattice spacing i.e. spins which are a large distance 
apart are uncorrelated. This is expected since in the 
thermodynamic limit, linir^oo pi^) = 0. The finite 
values of correlation function that appears in the figures 
is due to the finite size effect. 



In this model the simulation has been carried out for 
lattice sizes 5 x 5, 10 x 10, 15 x 15 and 20 x 20 using both 
WL and WLTM algorithms. Two dimensional random 
walk in two different spaces {E — S and E — p) was per- 
formed to determine the order parameter, susceptibility 
and correlation function apart from average energy, 
specific heat etc which can be determined using 1-d 
random walk. The minimum energy for all lattice sizes 
was 3. The upper limit of energies of these system sizes 



over which the simulation has been carried out were 50, 
200, 450 and 800. We deleted a small region at the lower 
end of energy to overcome the trapping of the random 
walker. The order parameter (correlation function) can 
have values between and 1 (—0.5 and +0.5). The 
whole energy and order parameter (correlation function) 
range is divided into a large number of bins of width 
de = 0.2 and d^ = 0.01 respectively. 
In figure 7, the density of states obtained from the 2-d 
random walk for 2-d XY model for different lattice sizes 
have been plotted against the energy per particle. We 
have used both WL and WLTM algorithms to obtain 
the data presented in figure 7. The order parameter 
{S{T)) for 2-d XY model of linear lattice sizes 5, 10, 
15 and 20 are plotted with temperature in figure 8. 
The system is known to possess no true long range 
order and a quasi-long-range-order disorder transition 
takes place due to unbinding of topological defects. 
Susceptibility of 2d XY model is also plotted as a 
function of temperature for four lattice sizes in figure 
9. It is observed that the height of the susceptibility 
peak increases with the increase of system sizes and also 
the position of susceptibility peak is shifted towards the 
lower temperature with the increase of system size. In 
figure 10 we have plotted the correlation function p(r, T) 
against temperature for ten values of r for the XY-model. 
These were obtained from the JDOS computed by using 
WL and WLTM methods. The same data is depicted 
in a different way in figure 10, where we have plotted 
p(r, T) against r for eight different values of temperature. 



VI. CONCLUSION 



We have presented the results of Monte Carlo simulation 
performed in the 1-d LL model and 2-d XY-model 
for the evaluation of joint density of states using 
the WL and WLTM algorithms. Agreement of the 
statistical averages of different quantities obtained by 
using the two algorithms is excellent. Calculation 
of JDOS for continuous spin models has been done 
earlier^. We have demonstrated in this work that, 
although computationally tedious, it is possible to 
use the WLTM method for the evaluation of the 
JDOS for a continuous spin model. This method may 
prove to be useful for future researchers who will need to 
generate the JDOS in a discrete or coninuous spin model. 
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